#ifndef GPUSolver_h
#define GPUSolver_h

// ============================================================
// Performance optimization: Disable ALL debug output
// Set to 0 to disable all printf/std::cout for maximum performance
// Set to 1 to enable debug output
// ============================================================
#define GPU_ENABLE_DEBUG_OUTPUT 0

// ============================================================
// Solver Runtime Information Output Control
// Set to 0 to disable runtime info (PCG iterations, timings, etc.)
// Set to 1 to enable solver runtime information
// Recommended: 0 for performance testing, 1 for debugging
// ============================================================
#define GPU_ENABLE_SOLVER_OUTPUT 0

// Solver output macro
#if GPU_ENABLE_SOLVER_OUTPUT
#define SOLVER_PRINTF(...) printf(__VA_ARGS__)
#else
#define SOLVER_PRINTF(...) ((void)0)
#endif

// ============================================================
// Mixed Precision Configuration
// Set to 1 to enable FP32 (single precision) on GPU side
// GPU uses FP32, CPU uses FP64, automatic conversion
// Expected speedup: 1.5-1.8x (theoretical 2x)
// ============================================================
#define GPU_USE_MIXED_PRECISION 1

#include <LinearSOE.h>
#include <Vector.h>
#include <vector>
#include <cuda_runtime.h>
#include <cusparse.h>
#include <cusolverSp.h>
#include <iostream>
#include "FastVectorOps.h"
#include "GPUElementSupport.h"

// Forward declarations
class Element;
class Domain;

// Forward declaration for GPUSolver
class GPUSolver;

// Debug output macro - defined outside class for global use
#if GPU_ENABLE_DEBUG_OUTPUT
#define GPU_DEBUG_PRINT(module, level, ...) \
	do { \
		if (GPUSolver::debugEnabled && \
			GPUSolver::moduleEnabled[static_cast<int>(module)] && \
			level <= GPUSolver::currentDebugLevel) { \
			std::cout << "[GPUSolver][" << GPUSolver::getModuleName(module) << "][" \
					  << GPUSolver::getLevelName(level) << "] " << __VA_ARGS__ << std::endl; \
		} \
	} while(0)
#else
#define GPU_DEBUG_PRINT(module, level, ...) ((void)0)
#endif

class GPUSolver : public LinearSOE {
public:
	GPUSolver();
	~GPUSolver();

	int getNumEqn(void) const override;
	int setSize(Graph& theGraph) override;
	int addA(const Matrix&, const ID&, double fact = 1.0) override;
	int addB(const Vector&, const ID&, double fact = 1.0) override;
	int setB(const Vector&, double fact = 1.0) override;
	void zeroA(void) override;
	void zeroB(void) override;
	void setX(int loc, double value) override;
	void setX(const Vector& x) override;
	const Vector& getX(void) override;
	const Vector& getB(void) override;
	double normRHS(void) override;
	int solve(void) override;

	// 非线性分析状态管理
	int commitState(void);  // 收敛后提交材料状态（trial -> committed）

	int sendSelf(int cTag, Channel& theChannel) override;
	int recvSelf(int cTag, Channel& theChannel, FEM_ObjectBroker& theBroker) override;

	void exportGlobalStiffnessMatrix(const char* filename);

	// ====== GPU Matrix Assembly Public Interface ======
	// GPU矩阵组装接口
	int formTangent_GPU();
	int formUnbalance_GPU();
	
	// 单元数据管理
	int initializeElementData();
	int updateElementData();

	// 分阶段初始化（修复DOF映射时序问题）
	int initializeGeometryAndMaterial();      // 阶段1：提取几何和材料数据（不依赖DOF_Group）
	int extractDOFMapping();                  // 阶段2：提取DOF映射（在formTangent_GPU首次调用时）
	int extractGeometryAndMaterialOnly();     // 辅助函数：仅提取几何和材料数据
	
	// GPU内存管理
	int allocateGPUMatrixMemory();
	void freeGPUMatrixMemory();
	
	// 单元支持检查
	bool isElementSupported(int elementClassTag);
	int getSupportedElementInfo();
	// ====== End of GPU Matrix Assembly Public Interface ======

	// ====== Debug Control Public Interface ======
	
	// Debug output levels
	enum class DebugLevel {
		ERROR = 0,
		WARNING = 1,
		INFO = 2,
		DEBUG = 3,
		VERBOSE = 4
	};
	
	// Debug output modules
	enum class DebugModule {
		MATRIX_ASSEMBLY = 0,
		GPU_COMPUTE = 1,
		MEMORY_MANAGEMENT = 2,
		DATA_TRANSFER = 3,
		ERROR_HANDLING = 4
	};
	
	// Global debug control
	static bool debugEnabled;
	
	// Module-specific debug control
	static bool moduleEnabled[5];  // 5 modules
	
	// Current debug level
	static DebugLevel currentDebugLevel;
	
	// Helper functions for debug output
	static const char* getModuleName(DebugModule module);
	static const char* getLevelName(DebugLevel level);
	
	static void enableDebug(bool enabled = true);
	static void disableDebug();
	static void setModuleDebugEnabled(int moduleIndex, bool enabled);
	static void setDebugVerbosity(int level); // 0=ERROR, 1=WARNING, 2=INFO, 3=DEBUG, 4=VERBOSE
	static void enableAllDebugModules();
	static void disableAllDebugModules();
	// ====== GPU Error Monitoring Public Interface ======
	// GPU错误统计和监控
	static int getGPUErrorCount();
	static int getCudaErrorCount();
	static int getCusparseErrorCount();
	static int getCusolverErrorCount();
	static void printGPUErrorStatistics();
	static void resetGPUErrorStatistics();
	// ====== End of GPU Error Monitoring Public Interface ======

	// ====== Performance Timing System Public Interface ======
	// 性能计时统计系统
	enum class TimerType {
		TOTAL_SOLVE = 0,
		MEMORY_TRANSFER = 1,
		GPU_COMPUTATION = 2,
		MATRIX_ASSEMBLY = 3,
		ERROR_HANDLING = 4,
		DATA_SYNCHRONIZATION = 5,
		ERROR_RECOVERY = 6
	};
	
	static void startTimer(TimerType timerType);
	static void stopTimer(TimerType timerType);
	static double getElapsedTime(TimerType timerType);
	static double getTotalElapsedTime(TimerType timerType);
	static void printTimingStatistics();
	static void resetTimingStatistics();
	static void enablePerformanceTiming(bool enabled = true);
	static void disablePerformanceTiming();
	static double getAverageTime(TimerType timerType);
	// ====== End of Performance Timing System Public Interface ======

	// ====== GPU Memory Monitoring Public Interface ======
	// GPU内存使用监控和优化系统
	static size_t getTotalGPUMemoryAllocated();
	static size_t getPeakGPUMemoryUsage();
	static int getGPUMemoryAllocationCount();
	static double getGPUMemoryEfficiency();
	static void printGPUMemoryStatistics();
	static void resetGPUMemoryStatistics();
	static void enableGPUMemoryMonitoring(bool enabled = true);
	static void disableGPUMemoryMonitoring();
	
	// GPU内存泄漏检测
	static bool detectGPUMemoryLeaks();
	static void reportGPUMemoryLeaks();
	
	// GPU内存状态查询
	static int getTotalAllocations();
	static size_t getCurrentAllocatedMemory();
	static size_t getPeakMemoryUsage();
	static void initializeGPUMemoryTracking();
	// ====== End of GPU Memory Monitoring Public Interface ======

	// ====== CPU-GPU Data Synchronization Public Interface ======
	// CPU-GPU数据同步系统
	int syncNodeDisplacementsToGPU(const Vector& displacements);
	int syncNodeVelocitiesToGPU(const Vector& velocities);  
	int syncNodeAccelerationsToGPU(const Vector& accelerations);
	int syncAllNodeDataToGPU(const Vector& displacements, 
	                         const Vector& velocities, 
	                         const Vector& accelerations);
	
	// 从GPU同步数据回CPU
	int syncNodeDisplacementsFromGPU(Vector& displacements);
	int syncAllNodeDataFromGPU(Vector& displacements, 
	                           Vector& velocities, 
	                           Vector& accelerations);
	
	// 数据同步状态查询
	bool isNodeDataSynchronized();
	void markNodeDataDirty();
	double getLastSyncTime();
	
	// 同步控制和优化
	void enableAutomaticSync(bool enabled = true);
	void disableAutomaticSync();
	bool isAutomaticSyncEnabled();
	// ====== End of CPU-GPU Data Synchronization Public Interface ======

	// ====== Material State Synchronization Public Interface ======
	// 材料状态变量同步系统
	int syncMaterialStatesToGPU();
	int syncMaterialStressToGPU(int elementTag, const double* stress, int stressSize);
	int syncMaterialStrainToGPU(int elementTag, const double* strain, int strainSize);
	int syncMaterialPropertiesToGPU(int elementTag, const double* properties, int propSize);
	
	// 批量材料状态同步
	int syncAllMaterialStatesToGPU();
	int syncMaterialStatesFromGPU();
	
	// 材料状态查询和控制
	bool isMaterialDataSynchronized();
	void markMaterialDataDirty();
	int getMaterialDataCount();
	double getMaterialSyncTime();
	
	// 材料状态缓冲区管理
	int initializeMaterialStateBuffers();
	// ====== End of Material State Synchronization Public Interface ======

	// ====== Element Geometry Update Public Interface ======
	// 单元几何/拓扑数据实时更新系统
	int updateElementGeometry();
	int updateElementCoordinates(int elementTag, const double* coords, int coordSize);
	int updateElementConnectivity(int elementTag, const int* connectivity, int nodeCount);
	int updateDeformedGeometry(const Vector& displacements);
	
	// 大变形支持
	int enableLargeDeformation(bool enabled = true);
	bool isLargeDeformationEnabled();
	int updateDeformationGradient();
	int computeUpdatedLagrangianGeometry();
	
	// 几何状态查询和控制
	bool isElementGeometryUpdated();
	void markElementGeometryDirty();
	double getGeometryUpdateTime();
	int getUpdatedElementCount();
	
	// 几何数据缓冲区管理
	int initializeGeometryBuffers();
	// ====== End of Element Geometry Update Public Interface ======

	// ====== Enhanced Debug Information Public Interface ======
	// 增强调试信息输出系统
	void printSynchronizationStatus();
	void printMemoryUsageSummary();
	void printPerformanceStatistics();
	void printDataTransferStatus();
	
	// 详细状态跟踪
	void logDataSyncEvent(const std::string& eventType, const std::string& details);
	void logMemoryOperation(const std::string& operation, size_t bytes, const std::string& bufferType);
	void logPerformanceEvent(const std::string& operation, double timeMs);
	void logGPUKernelLaunch(const std::string& kernelName, int elements, double timeMs);
	
	// 诊断和故障排除
	int validateAllBuffers();
	void printGPUDeviceInfo();
	void printSystemConfiguration();
	int generateDebugReport(const std::string& filename);
	
	// 调试输出控制
	void setDebugOutputFile(const std::string& filename);
	void enableDetailedLogging(bool enabled = true);
	bool isDetailedLoggingEnabled();
	
	// 增强调试系统集成功能
	void printCompleteSystemStatus();
	void enableComprehensiveDebugging();
	void disableComprehensiveDebugging();
	int runDiagnostics();
	void setupProductionDebugging();
	// ====== End of Enhanced Debug Information Public Interface ======

	// ====== GPU Vector Operations Public Interface ======
	// 高性能自定义GPU向量操作接口
	// 替换CUBLAS操作，实现8-18倍性能提升
	
	// 核心向量操作接口
	int fastVectorNorm2(int n, const double* d_x, double* result);        // 替换cublasDnrm2
	int fastVectorDot(int n, const double* d_x, const double* d_y, double* result); // 替换cublasDdot
	int fastVectorAxpy(int n, double alpha, const double* d_x, double* d_y);        // 替换cublasDaxpy
	int fastVectorCopy(int n, const double* d_x, double* d_y);             // 替换cublasDcopy
	int fastVectorScal(int n, double alpha, double* d_x);                  // 替换cublasDscal
	
	// CG求解器优化接口
	int optimizedCGVectorOps(int n, double* d_r, double* d_p, double* d_Ap, 
	                         double* d_x, double* scalars);
	int batchVectorOperations(int numOps, int* operations, void** params);
	
	// 向量操作管理器接口
	int initializeVectorOps(int maxVectorSize);
	int cleanupVectorOps();
	bool isVectorOpsInitialized() const;
	int setVectorOpOptimizationLevel(int level);
	
	// CUBLAS切换和性能对比接口
	int enableCUBLASFallback(bool enabled = true);
	int enableVectorOpPerformanceComparison(bool enabled = true);
	bool isCUBLASFallbackEnabled() const;
	void printVectorOpPerformanceComparison();
	double getVectorOpSpeedup(VectorOpTimerType operation);
	
	// 向量操作错误统计接口
	static int getVectorOpErrorCount();
	static int getVectorKernelLaunchErrorCount();
	static void printVectorOpErrorStatistics();
	static void resetVectorOpErrorStatistics();
	
	// 向量操作性能监控接口
	static void startVectorOpTimer(VectorOpTimerType timerType);
	static void stopVectorOpTimer(VectorOpTimerType timerType);
	static double getVectorOpElapsedTime(VectorOpTimerType timerType);
	static void printVectorOpTimingStatistics();
	static void resetVectorOpTimingStatistics();
	
	// 向量操作调试控制接口
	static void enableVectorOpDebug(VectorOpDebugModule module, bool enabled = true);
	static void setVectorOpDebugLevel(int level);
	static void enableAllVectorOpDebugModules();
	static void disableAllVectorOpDebugModules();
	// ====== End of GPU Vector Operations Public Interface ======

	// ====== GPU Assembly Control Public Interface ======
	// GPU组装控制接口
	bool isGPUAssemblyEnabled() const { return m_matrixOnGPU; }
	const char* getClassType() const override { return "GPUSolver"; }
	// ====== End of GPU Assembly Control Public Interface ======

	// ====== Preconditioner Control Public Interface ======
	// 预条件子控制接口
	void enablePreconditioner(bool enable);
	bool isPreconditionerEnabled() const;
	// ====== End of Preconditioner Control Public Interface ======

	// ====== Warm Start Control Public Interface ======
	// Warm Start控制接口
	void enableWarmStart(bool enable);
	bool isWarmStartEnabled() const;
	void resetWarmStart();  // 强制重置warm start数据
	void getWarmStartStatistics(int& warmStarts, int& coldStarts) const;
	// ====== End of Warm Start Control Public Interface ======

	// ====== End of Debug Control Public Interface ======

private:
	Vector X, B;
	std::vector<int> rowPtr, colInd; // CSR结构
	std::vector<double> val;         // CSR结构
	int size, nnz;
	int lastSize = 0; // 新增：上一次分配的size
	int lastNnz = 0;  // 新增：上一次分配的nnz

	// ===== Preconditioner相关成员 =====
	bool usePreconditioner;        // 是否使用预条件
	double* d_M_inv;               // Jacobi预条件矩阵（对角逆）GPU端
	double* d_z;                   // PCG中的z向量（z = M^(-1) * r）
	bool preconditionerBuilt;      // 预条件矩阵是否已构建

	// ===== Warm Start相关成员 =====
	bool warmStartEnabled;         // 是否启用warm start
	double* d_X_prev;              // 上一次的解（GPU端）
	bool hasWarmStartData;         // 是否有可用的warm start数据
	int warmStartCounter;          // warm start使用次数（用于统计）
	int coldStartCounter;          // cold start使用次数

	// ===== Mixed Precision (FP32) 相关成员 =====
#if GPU_USE_MIXED_PRECISION
	float* d_val_fp32;             // CSR矩阵值（FP32）
	float* d_B_fp32;               // 右端向量（FP32）
	float* d_X_fp32;               // 解向量（FP32）
	float* d_M_inv_fp32;           // 预条件矩阵（FP32）
	float* d_r_persistent_fp32;    // 残差向量（FP32）
	float* d_p_persistent_fp32;    // 搜索方向（FP32）
	float* d_Ap_persistent_fp32;   // A*p（FP32）
	float* d_z_fp32;               // 预条件向量（FP32）
	float* d_X_prev_fp32;          // Warm start解（FP32）
#endif

	// ===== 持久化GPU资源 =====
	static cublasHandle_t s_cublasHandle;
	static cusparseHandle_t s_cusparseHandle;
	static bool s_handlesInitialized;
	static int s_instanceCount;

	// 持久化工作向量（避免重复分配）
	double* d_r_persistent;
	double* d_p_persistent;
	double* d_Ap_persistent;
	void* dBuffer_persistent;
	size_t persistentSize;          // 当前分配的尺寸
	size_t persistentBufferSize;

	// 资源管理函数
	int initializePersistentResources();
	void freePersistentResources();

	// GPU相关
	int* d_rowPtr = nullptr;
	int* d_colInd = nullptr;
	double* d_val = nullptr;
	double* d_B = nullptr;
	double* d_X = nullptr;

	cusparseHandle_t handle = nullptr;
	cusparseMatDescr_t descrA = nullptr;
	cusolverSpHandle_t m_cusolverHandle = nullptr;  // 持久化的cuSolver句柄
	void* buffer = nullptr;
	size_t bufferSize = 0;

	// ====== 性能优化：跟踪GPU端数据状态（避免冗余传输）======
	bool m_matrixOnGPU = false;    // d_val是否已在GPU上组装完成
	bool m_rhsOnGPU = false;       // d_B是否已在GPU上组装完成
	// ====== End of GPU data state tracking ======

	void allocateGPU();
	void freeGPU();

	// ====== GPU Matrix Assembly Members ======
	// Single element data management
	int numElements;              // Total number of elements
	int maxDOFPerElement;         // Maximum DOF per element
	
	// GPU memory pointers for matrix assembly
	double* d_elementStiffness;   // All element stiffness matrix buffer
	double* d_elementResidual;    // All element residual vector buffer
	int* d_elementDOF;           // Element DOF mapping table
	int* d_elementTypes;         // Element type array
	double* d_elementData;       // Element geometry/material data

	// GPU memory for material states (Phase 2: Nonlinear materials)
	// 双状态缓冲机制：分离committed和trial状态以支持Newton迭代
	void* d_materialStates_committed;  // Committed material state (上次收敛的状态)
	void* d_materialStates_trial;      // Trial material state (当前试算状态)
	double* d_strainTotal;             // Committed total strain (上次收敛的应变)
	double* d_strainCurrent_trial;     // Trial current strain (当前试算应变)
	int numGaussPoints;                // Total number of integration points

	// GPU computation buffers
	double* d_tempMatrix;        // Temporary matrix computation space
	double* d_tempVector;        // Temporary vector computation space
	
	
	// Element statistics
	int numSupportedElements;    // Number of GPU-accelerated elements
	int numUnsupportedElements;  // Number of unsupported elements (CPU fallback)

	// CUDA streams for computation
	cudaStream_t computeStream;  // GPU computation stream
	cudaStream_t assemblyStream; // Matrix assembly stream

	// ====== Phase 1 Optimization: DOF Mapping Cache ======
	// DOF映射缓存（避免每步重新提取）
	bool m_dofMappingCached;           // DOF映射是否已缓存
	int m_numFE;                       // FE_Element数量
	int m_maxFEDOF;                    // FE_Element最大DOF数
	int* d_feDOFMapping;               // FE DOF映射（GPU缓冲区）
	int* d_feDOFCount;                 // 每个FE的DOF数量（GPU）
	int* d_feToElementMap;             // FE→Element索引映射（GPU）
	int extractAndCacheDOFMapping();   // 提取并缓存DOF映射（仅调用一次）
	// ====== End of Phase 1 Optimization ======

	// ====== Phase 2 Optimization: GPU Node Coordinate Update ======
	// GPU节点坐标更新（避免CPU端遍历单元）
	double* d_initialNodeCoords;       // 初始节点坐标（GPU）
	double* d_currentNodeCoords;       // 当前节点坐标（GPU）
	int m_numNodes;                    // 节点数量
	bool m_initialCoordsCached;        // 初始坐标是否已缓存
	int cacheInitialNodeCoordinates(); // 缓存初始节点坐标（仅调用一次）
	int updateNodeCoordinatesOnGPU();  // 在GPU上更新节点坐标

	// 持久化的节点位移GPU内存（避免频繁分配/释放）
	double* d_nodeDisplacementsCache = nullptr;
	size_t d_nodeDisplacementsCache_size = 0;
	// ====== End of Phase 2 Optimization ======

	// ====== CPU-GPU Data Synchronization Members ======
	// GPU节点数据缓冲区
	double* d_nodeDisplacements; // 节点位移向量 (GPU)
	double* d_nodeVelocities;    // 节点速度向量 (GPU)
	double* d_nodeAccelerations; // 节点加速度向量 (GPU)
	
	// 同步状态管理
	bool nodeDataSynchronized;   // 数据是否同步
	bool automaticSyncEnabled;   // 是否启用自动同步
	double lastSyncTime;         // 上次同步时间戳
	int nodeDataSize;            // 节点数据向量大小
	
	// 同步缓冲区
	double* h_tempNodeBuffer;    // CPU临时缓冲区，用于数据传输
	// ====== End of CPU-GPU Data Synchronization Members ======

	// ====== Material State Synchronization Members ======
	// GPU材料状态数据缓冲区
	double* d_materialStress;     // 材料应力状态GPU缓冲区
	double* d_materialStrain;     // 材料应变状态GPU缓冲区
	double* d_materialProperties; // 材料属性参数GPU缓冲区
	int* d_materialStateMap;      // 材料状态映射表GPU缓冲区
	
	// 材料状态管理
	bool materialDataSynchronized;  // 材料数据是否同步
	double materialSyncTime;        // 材料状态同步时间戳
	int maxMaterialStates;          // 最大材料状态数量
	int materialStateSize;          // 单个材料状态大小
	int materialPropertySize;       // 单个材料属性大小
	
	// 材料状态CPU缓冲区
	double* h_tempMaterialBuffer;   // CPU临时材料状态缓冲区
	std::vector<int> activeMaterials; // 活跃材料列表
	// ====== End of Material State Synchronization Members ======

	// ====== Element Geometry Update Members ======
	// GPU几何数据缓冲区
	double* d_elementCoordinates;    // 单元节点坐标GPU缓冲区
	double* d_deformedCoordinates;   // 变形后坐标GPU缓冲区
	double* d_deformationGradient;   // 变形梯度GPU缓冲区
	int* d_elementConnectivity;      // 单元连接关系GPU缓冲区
	double* d_referenceGeometry;     // 参考几何GPU缓冲区
	
	// 几何更新状态管理
	bool elementGeometryUpdated;     // 几何数据是否更新
	bool largeDeformationEnabled;    // 是否启用大变形模式
	double geometryUpdateTime;       // 几何更新时间戳
	int updatedElementCount;         // 已更新几何的单元数量
	int maxGeometryNodes;            // 最大几何节点数量
	
	// 几何数据CPU缓冲区
	double* h_tempGeometryBuffer;    // CPU临时几何缓冲区
	std::vector<int> updatedElements; // 已更新几何的单元列表
	// ====== End of Element Geometry Update Members ======

	// ====== Enhanced Debug Information Members ======
	// 增强调试信息状态
	bool detailedLoggingEnabled;      // 是否启用详细日志记录
	std::string debugOutputFile;     // 调试输出文件路径
	std::ofstream* debugFileStream;  // 调试文件流
	
	// 调试统计信息
	int totalDataSyncEvents;          // 总数据同步事件数
	int totalMemoryOperations;       // 总内存操作数
	int totalKernelLaunches;         // 总GPU核函数启动数
	double totalDataTransferTime;    // 总数据传输时间
	double totalKernelExecutionTime; // 总核函数执行时间
	
	// 详细跟踪缓冲区
	std::vector<std::string> syncEventLog;      // 同步事件日志
	std::vector<std::string> memoryOperationLog; // 内存操作日志
	std::vector<std::string> performanceEventLog; // 性能事件日志
	// ====== End of Enhanced Debug Information Members ======

	// ====== GPU Vector Operations Members ======
	// 高性能自定义GPU向量操作管理器
	VectorOpManager* vectorOpManager;    // 向量操作管理器实例
	
	// 向量操作配置和状态
	bool vectorOpsInitialized;           // 向量操作系统是否已初始化
	bool cublassFallbackEnabled;         // 是否启用CUBLAS回退
	bool vectorOpPerformanceComparison;  // 是否启用性能对比
	int vectorOpOptimizationLevel;       // 优化级别 (0=基础, 1=优化, 2=激进)
	int maxVectorSize;                   // 最大向量大小
	
	// 向量操作性能统计
	int totalVectorOperations;           // 总向量操作次数
	int fastVectorOperations;            // 使用自定义核函数的操作次数
	int cublasVectorOperations;          // 使用CUBLAS的操作次数
	double totalVectorOpTime;            // 总向量操作时间
	double totalCUBLASTime;              // 总CUBLAS操作时间
	
	// 向量操作错误统计
	int vectorOpErrorCount;              // 向量操作错误总数
	int vectorKernelErrors;              // 核函数启动错误
	int vectorMemoryErrors;              // 内存操作错误
	int vectorValidationErrors;          // 数值验证错误
	// ====== End of GPU Vector Operations Members ======

	// ====== GPU Matrix Assembly Private Methods ======
	// GPU核心计算函数
	int computeElementStiffness_GPU();
	int computeElementResidual_GPU();
	int assembleGlobalMatrix_GPU();
	int assembleGlobalVector_GPU();

	// ====== Preconditioner Private Methods ======
	int buildJacobiPreconditioner();
	// ====== End of Preconditioner Private Methods ======

	// 数据提取和准备
	int extractElementDataFromDomain();
	int extractMaterialProperties(Element* theElement, GPUElementData& elemData);

	// 数据传输管理
	int transferElementDataToGPU();
	int transferResultsFromGPU();
	
	// ====== CPU-GPU Data Synchronization Private Methods ======
	// 数据同步内部实现
	int allocateNodeDataBuffers();
	void freeNodeDataBuffers();
	int copyVectorToGPUBuffer(const Vector& source, double* d_target);
	int copyGPUBufferToVector(double* d_source, Vector& target);
	void updateSyncTimestamp();
	bool validateNodeDataSize(const Vector& data);
	// ====== End of CPU-GPU Data Synchronization Private Methods ======
	
	// ====== Material State Synchronization Private Methods ======
	// 材料状态同步内部实现
	int allocateMaterialStateBuffers();
	void freeMaterialStateBuffers();
	int copyMaterialDataToGPU(const double* source, double* d_target, int dataSize);
	int copyMaterialDataFromGPU(double* d_source, double* target, int dataSize);
	bool validateMaterialStateSize(int stateSize);
	void updateMaterialSyncTimestamp();
	int findMaterialSlot(int elementTag);
	int addActiveMaterial(int elementTag);
	void removeActiveMaterial(int elementTag);
	// ====== End of Material State Synchronization Private Methods ======

	// ====== Nonlinear Material State Management Private Methods ======
	// 非线性材料状态变量内存管理（用于GPU本构计算）
	int allocateNonlinearMaterialStates();
	void freeNonlinearMaterialStates();
	// ====== End of Nonlinear Material State Management Private Methods ======
	
	// ====== Element Geometry Update Private Methods ======
	// 几何更新内部实现
	int allocateGeometryBuffers();
	void freeGeometryBuffers();
	int copyGeometryDataToGPU(const double* source, double* d_target, int dataSize);
	int copyGeometryDataFromGPU(double* d_source, double* target, int dataSize);
	bool validateGeometrySize(int geometrySize);
	void updateGeometryTimestamp();
	int findElementGeometrySlot(int elementTag);
	int addUpdatedElement(int elementTag);
	void removeUpdatedElement(int elementTag);
	int computeDeformedCoordinates(const Vector& displacements);
	int updateElementJacobian(int elementTag);
	// ====== End of Element Geometry Update Private Methods ======
	
	// ====== Enhanced Debug Information Private Methods ======
	// 增强调试信息内部实现
	void initializeDebugLogging();
	void finalizeDebugLogging();
	void writeDebugHeader();
	void writeDebugFooter();
	std::string getCurrentTimestamp();
	std::string formatMemorySize(size_t bytes);
	std::string formatExecutionTime(double timeMs);
	void flushDebugOutput();
	void addToEventLog(std::vector<std::string>& log, const std::string& entry, int maxSize = 1000);
	// ====== End of Enhanced Debug Information Private Methods ======
	
	// ====== GPU Vector Operations Private Methods ======
	// 向量操作内部实现方法
	int initializeVectorOpManager();
	void cleanupVectorOpManager();
	int validateVectorOpParameters(int n, const double* d_x, const double* d_y = nullptr);
	
	// 向量操作路径选择
	int selectVectorOpPath(VectorOpTimerType opType, int vectorSize);
	int executeWithCUBLAS(VectorOpTimerType opType, void* params);
	int executeWithCustomKernel(VectorOpTimerType opType, void* params);
	
	// 性能监控和对比
	void recordVectorOpExecution(VectorOpTimerType opType, double executionTime, bool usedCUBLAS);
	void updateVectorOpStatistics(VectorOpTimerType opType, bool success, bool usedCUBLAS);
	double measureCUBLASOperation(VectorOpTimerType opType, void* params);
	double measureCustomKernelOperation(VectorOpTimerType opType, void* params);
	
	// 错误处理和验证
	void recordVectorOpError(int errorType, const std::string& operation, const std::string& details);
	bool validateVectorOpResults(VectorOpTimerType opType, const double* result, const double* expected, int n);
	int handleVectorOpFailure(VectorOpTimerType opType, int errorCode, bool fallbackToCUBLAS);
	
	// 调试和诊断
	void debugPrintVectorOp(VectorOpDebugModule module, int level, const std::string& message);
	std::string formatVectorOpStats(VectorOpTimerType opType);
	void logVectorOpEvent(const std::string& operation, int vectorSize, double time, bool usedCUBLAS);
	
	// 静态向量操作计时数据
	static double vectorOpStartTimes[6];     // 各操作开始时间戳
	static double vectorOpTotalTimes[6];     // 各操作累计时间
	static int vectorOpCounts[6];            // 各操作调用次数
	static bool vectorOpTimingEnabled;       // 向量操作计时是否启用
	
	// 静态向量操作调试控制
	static bool vectorOpDebugModules[5];     // 向量操作调试模块开关
	static int vectorOpDebugLevel;           // 向量操作调试级别
	
	// 静态CUBLAS性能对比数据
	static double lastCUBLASTime[6];         // 上次CUBLAS操作时间
	static double lastFastOpTime[6];         // 上次自定义操作时间
	static int operationComparisonCount[6];  // 各操作对比次数
	// ====== End of GPU Vector Operations Private Methods ======
	
	// ====== End of GPU Matrix Assembly Private Methods ======

};

#endif